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Abstract  Our  overall  goal  is  to  develop,  implement  and  transfer  accurate  new  nu¬ 
merical  methods  for  solving  moving  boundary  problems  in  materials  science.  We  have 
three  specific  objectives: 

•  Combine  semi-Lagrangian  time  stepping,  accurate  contouring  and  fast  geomet¬ 
ric  algorithms  to  develop  and  implement  accurate,  efficient  and  general  new 
methods  for  moving  sharp  interfaces. 

•  Develop  a  fast  modular  open  source  moving  interface  code  for  transfer  to  other 
researchers,  labs,  and  industry. 

•  Build  efficient,  accurate  and  general  integral  solvers  for  coupled  systems  of  par¬ 
tial  differential  equations  (PDEs)  modeling  common  material  phenomena,  and 
couple  these  solvers  to  our  modular  moving  interface  code. 

We  reached  a  milestone  with  respect  to  our  first  objective — an  accurate,  efficient 
and  general  moving  interface  method — with  the  semi-Lagrangian  contouring  method 
developed  in  References  [16-20].  We  are  now  finalizing  portable  C/C++  codes  for 
fast  adaptive  accurate  2D  contouring  and  general  2D  and  3D  moving  interfaces,  to 
attain  our  second  objective.  A  3D  code  has  been  developed,  tested  on  viscoelastic 
flow  [P0,P1,P2],  and  delivered  for  public  availability  on  SourceForge.  Our  2D  code 
and  a  preliminary  spectral  solver  based  on  Ewald  summation  has  been  implemented 
and  thoroughly  validated  on  2D  elastic  membranes  evolving  in  Stokes  incompressible 
flow  [P3] .  Work  is  ongoing  for  a  high-order  linearly-implicit  open-source  2D  code 
for  stiff  interfaces.  A  fast  and  extremely  robust  elliptic  solver  with  complex  inter¬ 
faces,  based  on  a  new  immersed  interface  approach  and  a  new  Krylov-accelerated 
multigrid  technique,  has  been  developed,  implemented  and  extensively  validated  on 
elliptic  problems  with  high-contrast  variable  coefficients  and  complex  interfaces  [P4]. 
Finally,  we  are  developing  fast  new  elliptic  PDE  solvers,  based  on  Ewald  summation 
[P5],  integral  equations  and  a  new  geometric  nonuniform  FFT,  to  complete  our  third 
objective. 
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Moving  interfaces  Our  work  on  moving  interface  problems  in  materials  science 
combines  fast  solvers  such  as  boundary  integral  methods,  deferred  correction  and 
Ewald  summation  with  fast  geometric  algorithms  and  semi-Lagrangian  implicit  rep¬ 
resentations  to  build  effective  new  numerical  methods. 

We  developed  an  implicit  boundary  integral  method  for  computing  periodic  den¬ 
drite  formation  in  the  symmetric  model  of  unstable  solidification  [7]  and  fast  algo¬ 
rithms  for  evaluating  heat  potentials  [2]  which  speeded  up  our  method  by  several 
orders  of  magnitude.  In  [6],  we  combined  the  boundary  integral  method  of  [7]  with 
fast  algorithms  from  [2,3,8]  and  the  level  set  method  of  [4]:  the  level  set  method  han¬ 
dled  topological  changes  effectively  while  fast  boundary  integral  techniques  ensured 
accuracy  and  efficiency  in  the  velocity  evaluation.  We  developed  and  analyzed  efficient 
and  accurate  new  vortex  methods  for  modeling  convection  in  the  melt  [5,13,14,15], 
together  with  new  error  analyses  [15]  and  quadrature  rules  [12]  for  general  integral 
equations.  In  the  interests  of  modeling  arbitrarily  complex  material  effects  in  the 
bulk,  PDE  solvers  based  on  Ewald  summation  and  analytic  preconditioning  were 
developed  in  [9,10,11].  Since  then,  we  have  focused  on  the  development  and  imple¬ 
mentation  of  highly  effective  new  numerical  methods  for  general  moving  interface 
problems  and  widely  applicable  related  computations  such  as  elliptic  PDE  solvers. 
We  summarize  three  projects  below:  fast  semi-Lagrangian  contouring  methods  for 
general  moving  interfaces  [16-20, P0-P2],  Ewald  summed  boundary  integral  methods 
for  Stokes  flow  with  moving  interfaces  [P3],  fast  robust  multigrid  solvers  for  elliptic 
problems  with  complex  interfaces  [P4],  and  fast  locally-corrected  spectral  methods 
for  elliptic  systems  in  periodic  geometry  [P5] . 

Moving  interface  problems  A  moving  interface  is  a  collection  T(f)  of  noninter¬ 
secting  oriented  closed  curves  (in  R2)  or  surfaces  (in  R3).  A  sufficiently  smooth 
moving  interface  has  an  outward  unit  normal  vector  N,  a  mean  curvature  C,  and  a 
normal  velocity  vector  VN.  A  moving  interface  problem  specifies  VN  as  a  functional 
of  r (t).  Examples  include  passive  transport  V  =  N  ■  F  where  F(x,t)  is  given,  geo¬ 
metric  motion  V  =  p{6)  —  q(6)C  where  cos  9  —  N  ■  x,  Ostwald  ripening  V  —  JA  , 
where  u  solves  the  Laplace  equation  A u  =  0  off  T(t)  and  u  —  C  on  T(t),  and  models 
for  crystal  growth  V  —  ^  ,  where  u  solves  the  heat  equation  ut  =  A u  off  T(f)  with 
the  geometric  boundary  condition  u  =  f(N ,  C,  V)  on  T(t). 


A  linearly-implicit  approach  Any  moving  interface  problem  can  be  reformulated 
as  a  PDE  for  a  function  tp  whose  zero  set  is  T(£).  The  normal,  curvature  and  velocity 


are  then 


V  (p 

"liWir 


C  =  —V  ■  N, 


VN 


<PtV<p 

IW 


Given  an  extension  of  DA  off  T(t)  to  a  globally  defined  function  W,  we  can  regard 
the  VN  formula  as  a  PDE  for  p: 


Pt  ~  W  ■  V<p  =  0. 
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We  solve  this  equation  on  an  adaptive  quadtree  mesh  to  eliminate  the  cost  of  going 
up  a  dimension.  Correct  viscosity  solutions  are  obtained  by  semi-Lagrangian  time 
stepping  with  exact  distancing  and  large  time  steps.  A  general  problem-independent 
velocity  extension  makes  our  method  modular  and  easy  to  apply. 

Semi-Lagrangian  methods  The  semi-Lagrangian  “CIR”  method  [1]  solves  tpt  — 
F(x,t)  ■  V (p  =  0  by  the  following  algorithm:  at  each  evaluation  point  x  at  the  new 
time  level  t  +  k,  move  x  back  with  velocity  F(x,  tn)  to  s  =  x  +  kF(x,  fn);  interpolate 
< p(x,tn )  to  the  point  s;  set  (p(x,tn+ 1)  equal  to  the  interpolated  value.  Our  second- 
order  time  stepping  scheme  couples  a  CIR  predictor  with  a  trapezoidal  corrector  using 
the  velocity  evaluated  from  the  CIR  approximation.  It  combines  the  unconditional 
stability  of  CIR  with  the  dramatically  reduced  dissipation  of  the  trapezoidal  rule. 
Interpolation  error  is  eliminated  by  exact  distance  finding  in  a  dynamic  quadtree 
data  structure.  The  tree  mesh  is  refined  with  a  functional  approach:  given  a  signed 
distance  function  <p(x,  tn),  we  build  a  tree  at  time  tn+1  —  tn+k  by  recursive  evaluation 
of  <p(x,tn+ 1)  =  < p(s,tn )  at  projected  points  s  =  x  +  kF(x,tn). 

Our  semi-Lagrangian  method  for  moving  interfaces  [19]  combines  efficient  exact 
quadtree-based  redistancing,  stable  second-order  semi-Lagrangian  time  stepping,  a 
modular  problem-independent  velocity  extension,  and  exact  ip  interpolation  in  the 
CIR  scheme.  The  velocity  extension  technique  evaluates  the  nearest-point  extension 
on  a  distance  tree,  builds  a  continuous  interpolant,  and  satisfies  a  maximum  principle. 
Our  method  resolves  and  moves  complex  interfaces  at  optimal  cost  with  time  steps 
unconstrained  by  numerical  stability.  It  is  a  “black-box”  method  for  moving  inter¬ 
faces,  which  accepts  the  interface  and  its  velocity  at  time  t  and  returns  the  evolved 
interface  one  time  step  later.  Such  methods  simplify  moving  interfaces,  because  the 
numerics  are  independent  of  the  physical  problem  driving  the  interfacial  motion.  Nu¬ 
merical  results  show  that  the  method  converges  to  correct  viscosity  solutions  even  for 
difficult  moving  interface  problems  involving  merging,  faceting,  transport,  nonlocality 
and  anisotropic  curvature-dependent  geometry. 

A  fast  semi-Lagrangian  contouring  method  [20]  General  moving  interface 
problems  are  solved  in  [20]  by  a  new  approach:  extract  the  moving  interface  from 
an  explicit  semi-Lagrangian  advection  formula  with  efficient  geometric  algorithms 
and  fast  accurate  contouring  techniques.  A  modular  adaptive  implementation  with 
fast  new  geometry  modules  computes  highly  accurate  solutions  to  moving  interface 
problems  involving  merging,  anisotropy,  faceting,  curvature,  dynamic  topology  and 
nonlocal  interactions  of  PDE  type.  Exact  geometric  algorithms  are  tuned  for  speed; 
velocity  evaluation  and  time  stepping  are  efficiently  decoupled  from  interface  reso¬ 
lution;  fast  new  contouring  techniques  dramatically  increase  overall  accuracy.  An 
efficient  adaptive  framework  combines  the  high  resolution  of  front  tracking  with  the 
topological  robustness  of  implicit  representations. 

Accurate  contouring  The  general  problem  of  finding  a  smooth  geometrically  con¬ 
strained  approximate  zero  set  of  a  function  which  can  be  evaluated  at  arbitrary  points 
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occurs  frequently  in  computational  science  and  requires  a  robust  general  contouring 
package.  An  ideal  contouring  package  would  accept  function  values  (and  derivatives 
if  available)  at  arbitrary  points  and  produce  a  piecewise-smooth  approximation  to 
the  zero  set  with  corners  where  necessary.  Geometric  constraints  such  as  bounds  on 
curvature  away  from  corners  are  vital  in  applications  such  as  computer-controlled 
machining,  and  pose  a  major  complication  for  existing  public-domain  contouring 
software.  The  PI,  with  a  graduate  student,  is  finalizing  new  open-source  C/C++ 
packages  for  constrained  piecewise-smooth  contouring  of  scattered  data  in  two  and 
three  dimensions. 

Ostwald  ripening  We  are  building  fast  nonlocal  velocity  evaluation  modules  for 
several  standard  moving  interface  problems  of  materials  science.  The  simplest  exam¬ 
ple  is  Ostwald  ripening,  which  models  the  growth  of  larger  solid  drops  by  evaporation 
from  smaller  drops  with  total  solid  volume  conserved.  The  velocity  V  is  the  normal 
derivative  JA  of  the  function  u  which  is  harmonic  off  T(f)  and  equal  to  the  cur¬ 
vature  on  T (t).  We  evaluate  this  nonlocal  velocity  by  solving  the  integral  equation 
of  classical  potential  theory  and  applying  the  Dirichlet  to  Neumann  operator.  The 
solution  u  is  a  double  layer  potential 

Dy{x)  =  X  d-mW^y)dy 

of  an  unknown  density  y  on  T  =  T(f),  with  K(x,y)  =  A  log  ||a;  —  y\\  the  free-space 
Green  function  for  the  two-dimensional  Laplace  equation.  The  density  y  solves  the 
integral  equation 

y(y)dy  =  C{x),  x  G  T. 

Once  y  is  found,  it  is  convenient  to  view  the  harmonic  function  u  as  the  real  part  of  an 
analytic  function  U .  The  Cauchy-Riemann  equations  then  yield  V  as  the  tangential 
derivative  of  the  imaginary  part  of  U,  which  is  easier  to  compute  than  the  normal 
derivative  of  the  real  part  u.  Discretization  of  this  formulation  is  highly  accurate  if 
the  interface  is  represented  by  equidistant  points  in  arclength.  Detailed  resolution 
of  the  interface  requires  many  points,  so  fast  algorithms  such  as  the  fast  multipole 
method  play  an  important  role. 

3D  viscoelastic  flow  with  interfaces  The  semi-Lagrangian  method  has  been 
implemented  in  3D  and  coupled  with  a  viscoelastic  fluid  simulator  to  produce  complex 
and  realistic  fluid  simulations  [P0,P1,P2],  A  particularly  useful  feature  of  the  semi- 
Lagrangian  approach  is  the  convenience  of  tracking  local  features  such  as  textures  on 
the  surface. 

Stokes  flow  with  interfaces  [P3]  Many  biological  moving  interface  problems 
(blood  flow,  cell  movement,  atherosclerosis)  involve  slow  viscous  flows  satisfying  the 
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incompressible  Stokes  equations 

—isAu  +  Vp  —  F,  V  •  u  —  0. 

Elastic  interfaces  produce  singular  forces  F  =  f8 r  where  <5r  is  a  measure  on  the 
interface  T.  We  have  developed  a  fast  solver  for  such  problems  by  combining  semi- 
Lagrangian  interface  evolution  with  a  fast  new  Ewald  summation  scheme  for  the 
Stokes  equations.  The  semi-Lagrangian  transport  of  interface  densities  allows  a 
straightforward  computation  of  stretching  energy,  while  the  new  fast  summation  tech¬ 
nique  unifies  several  well-known  local  correction  techniques  for  singular  integral  op¬ 
erators.  Our  fast  summation  technique  is  particularly  effective  for  Stokes  problems 
with  discontinuous  data  such  as  elastic  membrane  problems,  where  standard  methods 
encounter  great  difficulty.  Circularizing  and  oscillating  interfaces  under  Stokes  flow 
are  computed  accurately  to  one  part  per  thousand. 

Elliptic  solvers  with  complex  interfaces  [P4]  A  new  piecewise-polynomial  in¬ 
terface  method  (PIM)  for  discretizing  elliptic  problems  with  complex  interfaces  be¬ 
tween  high-contrast  materials  is  derived,  analyzed  and  tested.  A  Krylov-accelerated 
interface  multigrid  approach  (IMG)  solves  the  discretization  efficiently.  Stability  and 
convergence  are  proved  in  one  dimension,  while  an  extensive  array  of  numerical  ex¬ 
periments  with  complex  interfaces  and  large  coefficient  transitions  demonstrate  the 
accuracy,  efficiency  and  robustness  of  the  method  in  two  dimensions. 

Ewald  summation  for  elliptic  systems  [P5]  Our  Stokes  solver  generalizes  to 
solve  general  overdetermined  constant-coefficient  elliptic  systems 

AjUj  +  Aqu  =  f 

subject  to  periodic  boundary  conditions.  Thus  every  constant-coefficient  linear  ellip¬ 
tic  problem  can  be  solved  with  a  single  automatic  code,  eliminating  the  necessity  for 
many  specialized  codes.  Our  approach  yields  high-order  accuracy  even  for  solutions 
with  limited  smoothness,  and  combines  with  a  distributional  approach  to  arbitrary 
domains  to  yield  solvers  promising  arbitrary  order  accuracy  and  almost  optimal  effi¬ 
ciency. 
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